House Price - Modelagem
1 Objetivos
- Quais Modelos utilizar ?
- Modelos no Radar
- Regressão por Redes Bayesians
- Modelos Lineares Generalizados
- Modelos de Machine Learning
- Métricas de Performance
- RMSE
2 Conjunto de dados
model.train <- fread('../outputs/model.train.csv',
sep=",",
showProgress = FALSE)[,-1] %>%
data.frame(stringsAsFactors = T) %>%
select(Id,SalePrice,everything()) tipo <- lapply(model.train,class)
model.train[,unlist(tipo) != 'integer'] <- data.frame(apply(model.train[,unlist(tipo)!='integer'],2,factor))
model.train2.1 Importancia das variáveis
Obs: Esta parte do código ficará comentada para o markdown rodar mais rápido e para ilustrar o resultado irei apenas colocar a imagem da sugestão de seleção de variáveis feita pelo boruta, caso queiram rodar o código é só descomenta-lo.
selecao-de-variaveis
3 Ajustando um modelo aos dados
Como a distribuição da variável SalePrice é assimétrica com cauda a direita, podemos utilizar um MLG, modelos linerares generalizados, em particular a distribuição Gamma.
Em um primeiro momento utilizei as variáveis selecionadas pelo Boruta, 50 variáveis, porém após essa pré-seleção observei que algumas delas estavam com os p-valores muito alto então apliquei o método StepAic para continuar o processo de seleção o que me levou a configuração final abaixo.
Ainda é importante lembrar que dentro do contexto de negócio o processo de seleção de variáveis conta com o apoio dos especialistas da área, o que ajuda bastante, porém as vezes por questões intrínsecas do negócio uma variável que é rejeitada por algum método de seleção acaba permanecendo no modelo.
fit.model = glm(SalePrice ~ MSSubClass + LotArea + OverallQual +
OverallCond + YearBuilt + YearRemodAdd + MasVnrArea + BsmtUnfSF +
TotalBsmtSF + X2ndFlrSF + GrLivArea + BsmtFullBath + BedroomAbvGr +
KitchenAbvGr + TotRmsAbvGrd + Fireplaces + GarageYrBlt +
GarageCars + GarageArea + WoodDeckSF + MSZoning + Neighborhood +
HouseStyle + Exterior1st + MasVnrType + ExterQual + BldgType +
Foundation + BsmtQual + BsmtExposure + BsmtFinType1 + HeatingQC +
CentralAir + KitchenQual + Functional + SaleCondition,
data = model.train,family=Gamma(link=identity))
summary(fit.model)##
## Call:
## glm(formula = SalePrice ~ MSSubClass + LotArea + OverallQual +
## OverallCond + YearBuilt + YearRemodAdd + MasVnrArea + BsmtUnfSF +
## TotalBsmtSF + X2ndFlrSF + GrLivArea + BsmtFullBath + BedroomAbvGr +
## KitchenAbvGr + TotRmsAbvGrd + Fireplaces + GarageYrBlt +
## GarageCars + GarageArea + WoodDeckSF + MSZoning + Neighborhood +
## HouseStyle + Exterior1st + MasVnrType + ExterQual + BldgType +
## Foundation + BsmtQual + BsmtExposure + BsmtFinType1 + HeatingQC +
## CentralAir + KitchenQual + Functional + SaleCondition, family = Gamma(link = identity),
## data = model.train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.95166 -0.05437 0.00005 0.05192 0.40735
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.441e+05 1.166e+05 -7.238 8.15e-13 ***
## MSSubClass -1.158e+02 5.175e+01 -2.238 0.025386 *
## LotArea 1.066e+00 1.724e-01 6.186 8.48e-10 ***
## OverallQual 5.828e+03 7.163e+02 8.136 1.01e-15 ***
## OverallCond 4.953e+03 5.526e+02 8.963 < 2e-16 ***
## YearBuilt 2.999e+02 4.846e+01 6.187 8.41e-10 ***
## YearRemodAdd 6.302e+01 3.576e+01 1.762 0.078250 .
## MasVnrArea 8.576e+00 5.319e+00 1.612 0.107147
## BsmtUnfSF -8.395e+00 2.394e+00 -3.507 0.000471 ***
## TotalBsmtSF 2.694e+01 3.961e+00 6.801 1.64e-11 ***
## X2ndFlrSF 1.753e+01 5.209e+00 3.365 0.000791 ***
## GrLivArea 4.510e+01 4.026e+00 11.202 < 2e-16 ***
## BsmtFullBath 2.043e+03 1.281e+03 1.595 0.111087
## BedroomAbvGr -3.979e+03 9.753e+02 -4.080 4.79e-05 ***
## KitchenAbvGr -2.069e+04 3.885e+03 -5.327 1.19e-07 ***
## TotRmsAbvGrd 1.855e+03 7.252e+02 2.558 0.010640 *
## Fireplaces 2.953e+03 9.755e+02 3.027 0.002520 **
## GarageYrBlt 7.349e+01 3.656e+01 2.010 0.044637 *
## GarageCars 3.583e+03 1.578e+03 2.271 0.023344 *
## GarageArea 1.119e+01 5.874e+00 1.905 0.057027 .
## WoodDeckSF 1.118e+01 4.520e+00 2.474 0.013489 *
## MSZoningFV 4.486e+04 8.240e+03 5.445 6.30e-08 ***
## MSZoningRH 3.569e+04 6.307e+03 5.659 1.90e-08 ***
## MSZoningRL 3.813e+04 4.721e+03 8.077 1.61e-15 ***
## MSZoningRM 3.301e+04 4.107e+03 8.038 2.18e-15 ***
## NeighborhoodBlueste 1.952e+03 1.211e+04 0.161 0.871949
## NeighborhoodBrDale 1.163e+04 7.605e+03 1.529 0.126516
## NeighborhoodBrkSide 6.684e+03 7.113e+03 0.940 0.347590
## NeighborhoodClearCr -6.309e+03 7.791e+03 -0.810 0.418212
## NeighborhoodCollgCr -5.532e+03 5.999e+03 -0.922 0.356623
## NeighborhoodCrawfor 2.048e+04 6.862e+03 2.985 0.002893 **
## NeighborhoodEdwards -1.122e+04 6.348e+03 -1.767 0.077519 .
## NeighborhoodGilbert -6.596e+03 6.425e+03 -1.027 0.304795
## NeighborhoodIDOTRR 3.103e+03 7.630e+03 0.407 0.684285
## NeighborhoodMeadowV -4.235e+03 8.621e+03 -0.491 0.623340
## NeighborhoodMitchel -1.026e+04 6.615e+03 -1.551 0.121120
## NeighborhoodNAmes -9.638e+03 6.202e+03 -1.554 0.120452
## NeighborhoodNoRidge 3.150e+04 8.088e+03 3.895 0.000104 ***
## NeighborhoodNPkVill 1.112e+04 7.918e+03 1.405 0.160305
## NeighborhoodNridgHt 1.841e+04 6.706e+03 2.745 0.006148 **
## NeighborhoodNWAmes -1.001e+04 6.490e+03 -1.543 0.123171
## NeighborhoodOldTown -4.863e+03 7.151e+03 -0.680 0.496662
## NeighborhoodSawyer -6.549e+03 6.418e+03 -1.020 0.307786
## NeighborhoodSawyerW -3.635e+03 6.381e+03 -0.570 0.568993
## NeighborhoodSomerst 5.286e+03 8.398e+03 0.629 0.529185
## NeighborhoodStoneBr 3.403e+04 8.063e+03 4.220 2.63e-05 ***
## NeighborhoodSWISU 5.289e+03 7.194e+03 0.735 0.462392
## NeighborhoodTimber -7.486e+03 7.166e+03 -1.045 0.296382
## NeighborhoodVeenker 3.836e+03 9.242e+03 0.415 0.678159
## HouseStyle1.5Unf 1.125e+04 5.059e+03 2.223 0.026400 *
## HouseStyle1Story 8.747e+03 2.829e+03 3.092 0.002034 **
## HouseStyle2.5Unf 6.750e+03 5.962e+03 1.132 0.257769
## HouseStyle2Story -6.964e+02 2.345e+03 -0.297 0.766566
## HouseStyleSFoyer 1.105e+04 4.368e+03 2.529 0.011563 *
## HouseStyleSLvl 1.126e+04 3.737e+03 3.012 0.002646 **
## Exterior1stBrkComm -4.424e+04 9.799e+03 -4.515 6.96e-06 ***
## Exterior1stBrkFace 3.063e+03 4.264e+03 0.718 0.472777
## Exterior1stCemntBd -4.334e+03 5.919e+03 -0.732 0.464180
## Exterior1stHdBoard -1.072e+04 3.661e+03 -2.928 0.003477 **
## Exterior1stMetalSd -4.994e+03 3.448e+03 -1.448 0.147767
## Exterior1stPlywood -1.039e+04 4.110e+03 -2.528 0.011585 *
## Exterior1stStucco -3.908e+03 4.929e+03 -0.793 0.427954
## Exterior1stVinylSd -7.128e+03 3.542e+03 -2.013 0.044382 *
## Exterior1stWd Sdng -8.151e+03 3.516e+03 -2.318 0.020596 *
## Exterior1stWdShing -9.666e+03 4.692e+03 -2.060 0.039603 *
## MasVnrTypeBrkFace 8.424e+03 4.145e+03 2.032 0.042330 *
## MasVnrTypeNone 9.125e+03 4.195e+03 2.175 0.029802 *
## MasVnrTypeStone 1.269e+04 4.768e+03 2.661 0.007887 **
## ExterQualFa -3.218e+04 7.953e+03 -4.046 5.55e-05 ***
## ExterQualGd -2.008e+04 6.111e+03 -3.286 0.001046 **
## ExterQualTA -2.379e+04 6.284e+03 -3.785 0.000161 ***
## BldgType2fmCon 1.416e+04 7.529e+03 1.881 0.060277 .
## BldgTypeDuplex -6.343e+03 4.889e+03 -1.297 0.194725
## BldgTypeTwnhs -1.342e+04 6.609e+03 -2.030 0.042536 *
## BldgTypeTwnhsE -4.731e+03 5.954e+03 -0.795 0.427017
## FoundationCBlock 6.254e+03 1.972e+03 3.171 0.001556 **
## FoundationPConc 6.202e+03 2.196e+03 2.825 0.004808 **
## FoundationStone 1.810e+04 6.244e+03 2.899 0.003808 **
## FoundationWood -1.817e+04 1.108e+04 -1.640 0.101208
## BsmtQualFa -2.224e+04 4.616e+03 -4.818 1.64e-06 ***
## BsmtQualGd -2.166e+04 3.519e+03 -6.156 1.02e-09 ***
## BsmtQualTA -2.309e+04 3.890e+03 -5.937 3.79e-09 ***
## BsmtExposureGd 1.058e+04 2.595e+03 4.078 4.84e-05 ***
## BsmtExposureMn -3.661e+03 2.417e+03 -1.514 0.130172
## BsmtExposureNo -3.650e+03 1.846e+03 -1.978 0.048197 *
## BsmtFinType1BLQ 1.567e+02 1.796e+03 0.087 0.930457
## BsmtFinType1GLQ 3.790e+03 1.887e+03 2.009 0.044779 *
## BsmtFinType1LwQ -5.350e+03 2.323e+03 -2.303 0.021476 *
## BsmtFinType1Rec -1.618e+03 1.888e+03 -0.857 0.391794
## BsmtFinType1Unf -3.108e+03 2.016e+03 -1.542 0.123407
## HeatingQCFa 3.052e+03 3.008e+03 1.014 0.310573
## HeatingQCGd -1.430e+03 1.459e+03 -0.980 0.327152
## HeatingQCPo -2.351e+04 1.097e+04 -2.143 0.032317 *
## HeatingQCTA -2.675e+03 1.382e+03 -1.935 0.053188 .
## CentralAirY 3.736e+03 2.141e+03 1.745 0.081258 .
## KitchenQualFa -2.383e+04 4.824e+03 -4.940 8.93e-07 ***
## KitchenQualGd -2.266e+04 3.654e+03 -6.203 7.64e-10 ***
## KitchenQualTA -2.244e+04 3.756e+03 -5.975 3.03e-09 ***
## FunctionalMaj2 -4.983e+03 7.777e+03 -0.641 0.521827
## FunctionalMin1 5.310e+03 6.313e+03 0.841 0.400463
## FunctionalMin2 2.512e+03 6.218e+03 0.404 0.686356
## FunctionalMod 7.177e+02 8.199e+03 0.088 0.930256
## FunctionalTyp 1.384e+04 5.589e+03 2.476 0.013408 *
## SaleConditionAdjLand 2.315e+04 1.383e+04 1.674 0.094376 .
## SaleConditionAlloca -3.831e+03 6.013e+03 -0.637 0.524144
## SaleConditionFamily 2.856e+03 3.692e+03 0.773 0.439462
## SaleConditionNormal 3.334e+03 1.703e+03 1.958 0.050430 .
## SaleConditionPartial 1.120e+04 3.161e+03 3.543 0.000411 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.01073054)
##
## Null deviance: 189.616 on 1299 degrees of freedom
## Residual deviance: 13.611 on 1192 degrees of freedom
## AIC: 29343
##
## Number of Fisher Scoring iterations: 7
3.1 Teste Qui-Quadrado
round(1-pchisq(13.611,1192, ncp = 0, lower.tail = TRUE, log.p = FALSE),3)## [1] 1
O desvio do modelo foi de D(\(y\);\(\hat{\mu}\)) = 13.611, com 1192 graus de liberdade, que leva a P = 1,00 e indica um ajuste adequado.
4 Dados de teste
model.test <- fread('../outputs/model.test.csv',
sep=",",
showProgress = FALSE)[,-1] %>%
data.frame(stringsAsFactors = T) %>%
select(Id,everything()) tipo <- lapply(model.test,class)
model.test[,unlist(tipo) != 'integer'] <- data.frame(apply(model.test[,unlist(tipo)!='integer'],2,factor))
model.test4.1 Predição
predicao <- predict(fit.model,model.test)4.2 Plot Ecdf
Como não temos os valores da variável SalePrice no conjunto de dados de teste, podemos apenas verificar se a distribuição dos dados preditos esta de acordo com a distribuição dos dados reais da variável SalePrice do conjunto de dados de treino. A da distribuição impírica acumulada dos dados pode nos ajudar com isso.
latticeExtra::ecdfplot(~ fit.model$data$SalePrice + predicao,
auto.key=list(space='bottom',col = c('red','blue')),
col = c('red','blue'),
lwd = c(2,3),
xlab =" ",ylab = 'F(x)',
main = 'Distribuição Empírica Acumulada')Observe que há concordância entre os valores preditos e reais da variável SalePrice. Uma outra forma de verificar essa condordância é plotar as distribuições de densidade e frequência, como fiz abaixo.
predicao <- data.frame(predicao)p1 <- ggpubr::ggdensity(predicao, x = "predicao",
fill = "#0073C2FF", color = "black",
add = "mean", rug = TRUE) +
xlim(xlim = c(0,max(model.train$SalePrice)) ) +
ylim(ylim = c(0,8.0*10^(-6))) +
labs(title = 'Distribuição da densidade das predições') +
theme_dark()
p2 <- ggpubr::ggdensity(model.train, x = "SalePrice",
fill = "lightyellow", color = "black",
add = "mean", rug = TRUE) +
xlim(xlim = c(0,max(model.train$SalePrice)) ) +
ylim(ylim = c(0,8.0*10^(-6) )) +
labs(title = 'Distribuição da densidade da variável SalePrice') +
theme_dark()
gridExtra::grid.arrange(p1,p2,nrow = 1)Embora as distribuições estejam muito parecidas, os gráficos acima nos permite visualizar que os valores preditos assume uma variabilidade um pouco maior a partir do valor médio da distribuição.
p3 <- ggplot(predicao, aes(x = predicao)) +
geom_histogram(bins = 100,
color = "black",
fill = "#0073C2FF") +
xlim(xlim = c(0,max(model.train$SalePrice)) ) +
ylim(ylim = c(0,100)) +
labs(title = 'Distribuição de frequência das predições') +
theme_dark()
p4 <- ggplot(model.train, aes(x = SalePrice)) +
geom_histogram(bins = 100,
color = "black",
fill = "lightyellow") +
xlim(xlim = c(0,max(model.train$SalePrice)) ) +
ylim(ylim = c(0,100)) +
labs(title = 'Distribuição de frequência da variável SalePrice') +
theme_dark()
gridExtra::grid.arrange(p3,p4,nrow = 1)Se no conjunto de teste a variável SalePrice mantiver a mesma distribuição do conjunto de treino meu modelo terá boas chances de performar bem.
plot(fit.model$data$SalePrice,
fit.model$fitted.values,
col = c('blue','red'),
pch = 20,
type = 'p',
xlab = 'SalePrice',
ylab = 'Valores Ajustados',
main = 'SalesPrice vs Valores Preditos')
abline(0,1)
legend('topleft',legend = c('SalePrice','Valores ajustados'),col = c('blue','red'), lty = c(2,2) )Este plot mostra que nosso modelo é eficaz para valores até \(4e+05\), isso ocorre por existirem poucas casas com valores muito alto no conjunto de dados, pois a variabilidade tem convergência assintótica, ou seja, \(\sigma_{_{n \implies \infty}} \implies 0\). Uma alternativa visando melhorar a performance do modelo para valores acima de \(4e+05\) é criar variáveis que ajudam a discriminar a variável resposta neste segmento. Através de uma clusterização acho que isso seria possível.
4.3 Resíduos
plot(fit.model$residuals, pch = 20, type = 'p', main = 'Resíduos')Os resíduos estão bem distribuídos em torno de zero e isso é muito bom! Considerando que o ajuste está bom, é interessante notar que os resíduos que estão muito acima de zero podem ser casos em que o valor do imóvel esta fora do praticado pelo mercado e isto implica que dificilmente um imóvel nesta situação seria facilmente vendido fazendo com que a equipe de vendas perdesse tempo e energia. Já os imóveis com preços abaixo de zero seriam facilmente vendidos porém o lucro da venda seria baixo.
5 Rsquare
rss <- sum((fit.model$residuals)^2) ## soma dos quadrados dos resíduos
tss <- sum((fit.model$data$SalePrice - mean(fit.model$data$SalePrice))^2) ## soma total dos quadrados
(r.square <- 1-rss/tss)## [1] 0.8942537
Este é um bom resultado para uma primeira iteração, o modelo consegue explicar cerca de 89% da variabilidade dos dados.
6 RMSE
sqrt(mean((fit.model$residuals)^2))## [1] 25378.23
Com esse RMSE somente o kaggle me dirá se está bom ou não.
7 Gráficos de Dispersão
layout(matrix(c(1,1,2,2,3,3), ncol = 2, nrow = 3, byrow = T))
plot(fit.model$data$SalePrice,
col = 'black',
pch = 20,
type = 'p',
ylim = c(0,7*10^5),
xlab = 'SalePrice',
ylab = '',
main = 'Valores reais de SalesPrice')
plot(fit.model$fitted.values,
col = 'blue',
pch = 20,
type = 'p',
ylim = c(0,7*10^5),
xlab = 'SalePrice',
ylab = '',
main = 'Valores Ajustados de SalesPrice')
plot(predicao$predicao,
col = 'red',
pch = 20,
type = 'p',
ylim = c(0,7*10^5),
xlab = 'SalePrice',
ylab = '',
main = 'Valores Preditos com df.test para SalesPrice')8 QQ-Plot
O QQ-plot pode nos ajudar a ter uma idéia de como será o resultados da predição, aparemtemente será bom.
qqplot(sort(fit.model$data$SalePrice)[1:1318],
sort(predicao$predicao),
main = 'QQ-plot',
ylab = 'Predição',
xlab = 'SalePrice',
col = c('blue','red') )
abline(0,1)9 Exportando para submissão no Kaggle
submit.glm <- data.frame(Id = model.test$Id,SalePrice = predicao$predicao, stringsAsFactors = F)
write.csv(submit.glm, file="../outputs/saleprice_glm.csv",row.names = FALSE,quote=FALSE)10 Conclusão
Neste primeiro momento os resultados são razoáveis.